#R version 4.2.1 (2022-06-23)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 10 x64 (build 19044)

library(dplyr) #version 1.0.9
library(readr) #version 2.1.2
library(stringr) #version 1.4.0
library(zoo) #version 1.8-10
library(sandwich) #version 3.0-2
library(lmtest) #version 0.9-40
library(prediction) #version 0.3.14
library(ggpubr) #version 0.4.0
library(estimatr) #version 1.0.0
library(margins) #version 0.3.26
library(broom) #version 1.0.0
library(scales) #version 1.2.0


`%notin%` <- function(x,y) !(x %in% y) 

#GWF Models - Table 1####
d1 <- read_csv("data/gwf_fa_pi.csv") 
d2 <- read_csv("data/gwf_irt.csv") %>% 
  rename(gwf_casename = case) %>%
  select(gwf_casename, transition, gwf_irt) 

gwf <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("gwf_casename", "year")) %>% 
  left_join(d2, by = c("gwf_casename")) %>% 
  group_by(gwf_casename) %>% 
  mutate_at(.vars = vars(gwf_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, gwf_casename, gwf_regimetype, v2xps_party, v2pssunpar, gwf_incumbent_pi, gwf_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, gwf_priorduration, coups, gwf_number, EFindex, theta_mean) %>%
  mutate(gwf_incumbent_pi = ifelse(is.na(gwf_incumbent_pi), v2xps_party, gwf_incumbent_pi)) %>% 
  rename(incumbent_pi = gwf_incumbent_pi) %>% 
  filter(transition == year, gwf_casename %notin% c("Bolivia 43-46", "Thailand 57-73", "Pakistan 58-71", "Haiti 88-90",
                                                    "Sudan 58-64", "Ecuador 63-66", "Venezuela 48-58")) %>% 
  ungroup() %>% 
  unique()

m1 <- lm(gwf_irt ~ incumbent_pi, data = gwf) 
se1    <- sqrt(diag(vcovHC(m1, type = "HC1")))

m2 <- lm(gwf_irt ~ incumbent_pi + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf) 
se2    <- sqrt(diag(vcovHC(m2, type = "HC1")))

m3 <- lm(gwf_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf) 
se3    <- sqrt(diag(vcovHC(m3, type = "HC1")))

#strength
d1 <- read_csv("data/gwf_fa_strength.csv") 
d2 <- read_csv("data/gwf_irt.csv") %>% 
  rename(gwf_casename = case) %>%
  dplyr::select(gwf_casename, transition, gwf_irt) 

gwf <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("gwf_casename", "year")) %>% 
  left_join(d2, by = c("gwf_casename")) %>% 
  group_by(gwf_casename) %>% 
  mutate_at(.vars = vars(gwf_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

gwf$v2pssunpar <- scales::rescale(gwf$v2pssunpar, to = c(0, 1), from = range(gwf$v2pssunpar, na.rm = TRUE, finite = TRUE))

gwf %<>%
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, gwf_casename, gwf_regimetype, v2xps_party, v2pssunpar, gwf_incumbent_strength, 
                gwf_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, gwf_priorduration, coups, gwf_number, EFindex, theta_mean) %>%
  mutate(gwf_incumbent_strength = ifelse(is.na(gwf_incumbent_strength), v2pssunpar, gwf_incumbent_strength)) %>% 
  rename(incumbent_strength = gwf_incumbent_strength) %>% 
  #filter outliers
  filter(transition == year, gwf_casename %notin% c("Bolivia 43-46", "Thailand 57-73", "Pakistan 58-71", "Haiti 88-90",
                                                    "Sudan 58-64", "Ecuador 63-66", "Venezuela 48-58")) %>% 
  ungroup() %>% 
  unique() 

m4 <- lm(gwf_irt ~ incumbent_strength, data = gwf) 
se4    <- sqrt(diag(vcovHC(m4, type = "HC1")))

m5 <- lm(gwf_irt ~ incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf) 
se5    <- sqrt(diag(vcovHC(m5, type = "HC1")))

m6 <- lm(gwf_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf) 
se6    <- sqrt(diag(vcovHC(m6, type = "HC1")))


#produce table
stargazer(m1, m2, m3, m4, m5, m6,  se = list(se1, se2, se3, se4, se5, se6), omit.stat=c("f", "ser"), omit = "Constant", 
          order = c(1, 17), covariate.labels = c("Institutionalization", "Strength",  "GINI", "Rural Inequality", 
                                                 "Military Spending", "Area", "Urban Pop", "GDP PC", 
                                                 "Region", "Resource Wealth", "Islam", "Education", "Duration", "Coups", "Regimes", "Ethnic Frac", "Repression"),
          column.labels   = c("Institutionalization", "Strength"), 
          column.separate = c(3, 3))


#PAR Models - Table 2####
#institutionalization
d1 <- read_csv("data/svolik_fa_pi.csv") 
d2 <- read_csv("data/svolik_irt.csv") %>% 
  rename(svolik_case = case) %>%
  select(svolik_case, transition, svolik_irt)

svolik <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("svolik_case", "year")) %>% 
  left_join(d2, by = c("svolik_case")) %>% 
  group_by(svolik_case) %>% 
  mutate_at(.vars = vars(svolik_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, svolik_case, svolik_regime, v2xps_party, v2pssunpar, svolik_incumbent_pi, svolik_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea,
         svolik_priorduration, coups, svolik_number, EFindex, theta_mean) %>%
  mutate(svolik_incumbent_pi = ifelse(is.na(svolik_incumbent_pi), v2xps_party, svolik_incumbent_pi)) %>% 
  rename(incumbent_pi = svolik_incumbent_pi)  %>% 
  filter(transition == year, svolik_case %notin% c("Bolivia 1946-1946", "Thailand 1971-1972", "Pakistan 1958-1972", "Haiti 1986-1990", 
                                                   "Sudan 1959-1964", "Ecuador 1963-1965", "Haiti 1950-1956",
                                                   "Comoros 1999-2001", "Venezuela 1948-1949", "Burkina Faso (Upper Volta) 1966-1971", "Peru 1949-1956"))%>% 
  ungroup() %>% 
  unique()

m1 <- lm(svolik_irt ~ incumbent_pi, data = svolik)
se1 <- sqrt(diag(vcovHC(m1, type = "HC1")))

m2 <- lm(svolik_irt ~ incumbent_pi + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik)
se2    <- sqrt(diag(vcovHC(m2, type = "HC1")))

m3 <- lm(svolik_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik)
se3    <- sqrt(diag(vcovHC(m3, type = "HC1")))

#strength
d1 <- read_csv("data/svolik_fa_strength.csv") 
d2 <- read_csv("data/svolik_irt.csv") %>% 
  rename(svolik_case = case) %>%
  dplyr::select(svolik_case, transition, svolik_irt)

svolik <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("svolik_case", "year")) %>% 
  left_join(d2, by = c("svolik_case")) %>% 
  group_by(svolik_case) %>% 
  mutate_at(.vars = vars(svolik_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

svolik$v2pssunpar <- scales::rescale(svolik$v2pssunpar, to = c(0, 1), from = range(svolik$v2pssunpar, na.rm = TRUE, finite = TRUE))

svolik %<>%  
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar))%>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, svolik_case, svolik_regime, v2xps_party, v2pssunpar, svolik_incumbent_strength, 
                svolik_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea,
                svolik_priorduration, coups, svolik_number, EFindex, theta_mean) %>%
  mutate(svolik_incumbent_strength = ifelse(is.na(svolik_incumbent_strength), v2pssunpar, svolik_incumbent_strength)) %>% 
  rename(incumbent_strength = svolik_incumbent_strength) %>% 
  filter(transition == year, svolik_case %notin% c("Bolivia 1946-1946", "Thailand 1971-1972", "Pakistan 1958-1972", "Haiti 1986-1990", 
                                                   "Sudan 1959-1964", "Ecuador 1963-1965", "Haiti 1950-1956",
                                                   "Comoros 1999-2001", "Venezuela 1948-1949", "Burkina Faso (Upper Volta) 1966-1971", "Peru 1949-1956"))%>% 
  ungroup() %>% 
  unique()

m4 <- lm(svolik_irt ~ incumbent_strength, data = svolik)
se4 <- sqrt(diag(vcovHC(m4, type = "HC1")))

m5 <- lm(svolik_irt ~ incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik)
se5    <- sqrt(diag(vcovHC(m5, type = "HC1")))

m6 <- lm(svolik_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik)
se6    <- sqrt(diag(vcovHC(m6, type = "HC1")))

stargazer(m1, m2, m3, m4, m5, m6,  se = list(se1, se2, se3, se4, se5, se6), omit.stat=c("f", "ser"), omit = "Constant", 
          order = c(1, 17), covariate.labels = c("Institutionalization", "Strength",  "GINI", "Rural Inequality", 
                                                 "Military Spending", "Area", "Urban Pop", "GDP PC", 
                                                 "Region", "Resource Wealth", "Islam", "Education", "Duration", "Coups", "Regimes", "Ethnic Frac", "Repression"),
          column.labels   = c("Institutionalization", "Strength"), 
          column.separate = c(3, 3))




#WTH Models - Table 3####
#institutionalization##
d1 <- read_csv("data/wth_fa_pi.csv") 
d2 <- read_csv("data/wth_irt.csv") %>% 
  rename(wth_case = case) %>%
  select(wth_case, transition, wth_irt)

wth <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("wth_case", "year")) %>% 
  left_join(d2, by = c("wth_case")) %>% 
  group_by(wth_case) %>% 
  mutate_at(.vars = vars(wth_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, wth_case, wth_regime, v2xps_party, v2pssunpar, wth_incumbent_pi, wth_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
         wth_priorduration, coups, wth_number, EFindex, theta_mean) %>%
  mutate(wth_incumbent_pi = ifelse(is.na(wth_incumbent_pi), v2xps_party, wth_incumbent_pi)) %>% 
  rename(incumbent_pi = wth_incumbent_pi)  %>% 
  filter(transition == year, wth_case %notin% c("Haiti 1986-1989", "Comoros 1999-2001")) %>% 
  ungroup() %>% 
  unique()

m1 <- lm(wth_irt ~ incumbent_pi, data = wth) 
se1    <- sqrt(diag(vcovHC(m1, type = "HC1")))

m2 <- lm(wth_irt ~ incumbent_pi + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth)  
se2    <- sqrt(diag(vcovHC(m2, type = "HC1")))

m3 <- lm(wth_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth) 
se3    <- sqrt(diag(vcovHC(m3, type = "HC1")))

#strength models #
d1 <- read_csv("data/wth_fa_strength.csv") 
d2 <- read_csv("data/wth_irt.csv") %>% 
  rename(wth_case = case) %>%
  dplyr::select(wth_case, transition, wth_irt)

wth <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("wth_case", "year")) %>% 
  left_join(d2, by = c("wth_case")) %>% 
  group_by(wth_case) %>% 
  mutate_at(.vars = vars(wth_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

wth$v2pssunpar <- scales::rescale(wth$v2pssunpar, to = c(0, 1), from = range(wth$v2pssunpar, na.rm = TRUE, finite = TRUE))

wth %<>%
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, wth_case, wth_regime, v2xps_party, v2pssunpar, wth_incumbent_strength, 
                wth_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
                wth_priorduration, coups, wth_number, EFindex, theta_mean) %>%
  mutate(wth_incumbent_strength = ifelse(is.na(wth_incumbent_strength), v2pssunpar, wth_incumbent_strength)) %>% 
  rename(incumbent_strength = wth_incumbent_strength) %>% 
  filter(transition == year, wth_case %notin% c("Haiti 1986-1989", "Comoros 1999-2001")) %>% 
  ungroup() %>% 
  unique() 

m4 <- lm(wth_irt ~ incumbent_strength, data = wth) 
se4    <- sqrt(diag(vcovHC(m4, type = "HC1")))

m5 <- lm(wth_irt ~ incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth) 
se5    <- sqrt(diag(vcovHC(m5, type = "HC1")))

m6 <- lm(wth_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth)
se6    <- sqrt(diag(vcovHC(m6, type = "HC1")))

stargazer(m1, m2, m3, m4, m5, m6,  se = list(se1, se2, se3, se4, se5, se6), omit.stat=c("f", "ser"), omit = "Constant", 
          order = c(1, 17), covariate.labels = c("Institutionalization", "Strength",  "GINI", "Rural Inequality", 
                                                 "Military Spending", "Area", "Urban Pop", "GDP PC", 
                                                 "Region", "Resource Wealth", "Islam", "Education", "Duration", "Coups", "Regimes", "Ethnic Frac", "Repression"),
          column.labels   = c("Institutionalization", "Strength"), 
          column.separate = c(3, 3))
#DD Models - Table 4####
#institutionalization
d1 <- read_csv("data/dd_fa_pi.csv") 
d2 <- read_csv("data/dd_irt.csv") %>% 
  rename(dd_case = case) %>%
  select(dd_case, transition, dd_irt)

dd <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("dd_case", "year")) %>% 
  left_join(d2, by = c("dd_case")) %>% 
  group_by(dd_case) %>% 
  mutate_at(.vars = vars(dd_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, dd_case, dd_regime, v2xps_party, v2pssunpar, dd_incumbent_pi, dd_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
         dd_priorduration, coups, dd_number, EFindex, theta_mean) %>%
  mutate(dd_incumbent_pi = ifelse(is.na(dd_incumbent_pi), v2xps_party, dd_incumbent_pi)) %>% 
  rename(incumbent_pi = dd_incumbent_pi)  %>% 
  filter(transition == year, dd_case %notin% c("Bolivia 1946-1946", "Thailand 1946-1972", "Pakistan 1958-1970", "Haiti 1986-1989", 
                                               "Sudan 1958-1963", "Ecuador 1963-1965", "Haiti 1950-1955",
                                               "Comoros 1999-2003", "Venezuela 1948-1958")) %>% 
  ungroup() %>% 
  unique()

m1 <- lm(dd_irt ~ incumbent_pi, data = dd)
se1 <- sqrt(diag(vcovHC(m1, type = "HC1")))

m2 <- lm(dd_irt ~ incumbent_pi + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd)
se2    <- sqrt(diag(vcovHC(m2, type = "HC1")))

m3 <- lm(dd_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd)
se3    <- sqrt(diag(vcovHC(m3, type = "HC1")))


#strength
d1 <- read_csv("data/dd_fa_strength.csv") 
d2 <- read_csv("data/dd_irt.csv") %>% 
  rename(dd_case = case) %>%
  dplyr::select(dd_case, transition, dd_irt)

dd <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("dd_case", "year")) %>% 
  left_join(d2, by = c("dd_case")) %>% 
  group_by(dd_case) %>% 
  mutate_at(.vars = vars(dd_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

dd$v2pssunpar <- scales::rescale(dd$v2pssunpar, to = c(0, 1), from = range(dd$v2pssunpar, na.rm = TRUE, finite = TRUE))

dd %<>%
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, dd_case, dd_regime, v2xps_party, v2pssunpar, dd_incumbent_strength, 
                dd_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
                dd_priorduration, coups, dd_number, EFindex, theta_mean) %>%
  mutate(dd_incumbent_strength = ifelse(is.na(dd_incumbent_strength), v2pssunpar, dd_incumbent_strength)) %>% 
  rename(incumbent_strength = dd_incumbent_strength) %>% 
  filter(transition == year, dd_case %notin% c("Bolivia 1946-1946", "Thailand 1946-1972", "Pakistan 1958-1970", "Haiti 1986-1989", 
                                               "Sudan 1958-1963", "Ecuador 1963-1965", "Haiti 1950-1955",
                                               "Comoros 1999-2003", "Venezuela 1948-1958")) %>% 
  ungroup() %>% 
  unique() 

m4 <- lm(dd_irt ~ incumbent_strength, data = dd)
se4 <- sqrt(diag(vcovHC(m4, type = "HC1")))

m5 <- lm(dd_irt ~ incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd)
se5    <- sqrt(diag(vcovHC(m5, type = "HC1")))

m6 <- lm(dd_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd)
se6    <- sqrt(diag(vcovHC(m6, type = "HC1")))

stargazer(m1, m2, m3, m4, m5, m6,  se = list(se1, se2, se3, se4, se5, se6), omit.stat=c("f", "ser"), omit = "Constant", 
          order = c(1, 17), covariate.labels = c("Institutionalization", "Strength",  "GINI", "Rural Inequality", 
                                                 "Military Spending", "Area", "Urban Pop", "GDP PC", 
                                                 "Region", "Resource Wealth", "Islam", "Education", "Duration", "Coups", "Regimes", "Ethnic Frac", "Repression"),
          column.labels   = c("Institutionalization", "Strength"), 
          column.separate = c(3, 3))


